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We present a brief survey of methods that utihze computer simulations and quantum and statis- 
tical mechanics in the analysis of electrochemical systems. The methods, Molecular Dynamics and 
Monte Carlo simulations and quantum-mechanical density-functional theory, are illustrated with 
examples from simulations of lithium-battery charging and electrochemical adsorption of bromine 
on single-crystal silver electrodes. 



I. INTRODUCTION 

The interface between a solid electrode and a liquid electrolyte is a complicated many-particle system, in which the 
electrode ions and electrons interact with solute ions and solvent ions or molecules through several channels of inter- 
action, including forces due to quantum-mechanical exchange, electrostatics, hydrodynamics, and elastic deformation 
of the substrate. Over the last few decades, surface electrochemistry has been revolutionized by new techniques that 
enable atomic-scale observation and manipulation of solid-liquid interfaces [11, , yielding novel methods for materials 
analysis, synthesis, and modification. This development has been paralleled by equally revolutionary developments 
in computer hardware and algorithms that by now enable simulations with millions of individual particles 3] , so that 
there is now significant overlap between system sizes that can be treated computationally and experimentally. 

In this chapter, we discuss some of the methods available to study the structure and dynamics of electrode- 
electrolyte interfaces using computers and techniques based on quantum and statistical mechanics. These methods 
are illustrated by some recent applications. The rest of the chapter is organized as follows. In Sec. [ITl we present 
fully three-dimensional, continuum simulations by Molecular Dynamics (MD) of ion intercalation during charging of 
Lithium- ion batteries. In Sec. lIIIl we discuss the simplifications that are possible by mapping a chemisorption problem 
onto an effective lattice-gas (LG) Hamiltonian , and in Sec. llVl we demonstrate how input parameters for a statistical- 
mechanical LG model can be estimated from quantum-mechanical density-functional theory (DFT) calculations. 
Section |V] is devoted to a discussion of Monte Carlo (MC) simulations, both for equilibrium problems (Sec. IV A]) and 
for dynamics (Sec. IV B]) . As an example of the latter, we present in Sec. lVIl a simulational demonstration of a method 
to classify surface phase transitions in adsorbate systems, which is an extension of standard cyclic voltammetry (CV): 
the Electrochemical First-order Reversal Curve (EC-FORC) method. A concluding summary is given in Sec. IVIII 
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FIG. 1: (a) Snapshot of the model system containing four graphite sheets, two PFg ions and ten Li"*" ions (spheres), solvated in 
69 propylene carbonate and 87 ethylene carbonate molecules after reaching constant volume in the NPT ensemble, (b) Snapshot 
after 200 ns MD simulation. The ensemble is the NVT ensemble. The system has periodic boundary conditions and is simulated 
at one atm and 300 K. Top view, perpendicular to the plane of the graphite sheets. 

II. MOLECULAR DYNAMICS SIMULATIONS OF ION INTERCALATION IN LITHIUM BATTERIES 

The charging process in Lithium-ion batteries is marked by the intercalation of Lithium ions into the graphite anode 
material. Here we present MD simulations of this process and suggest a new charging method that has the potential 
for shorter charging times, as well as the possibility of providing higher power densities. 

A. Molecular Dynamics and Model System 

Molecular Dynamics is based on solving the classical equations of motion for a system of N atoms interacting 
through forces derived from a potential-energy function 0, H, H, 0, Hi • From the potential energy , the force on the 
ith atom, Fi, is calculated. Thus, the equation of motion is 



dEp 



' dt 



(1) 



where rj, Vi, and mi are the position, velocity, and mass of the ith atom, respectively. Consequently, the quality of the 
simulations strongly depends on the ability of the classical force field to reasonably describe the atomistic behavior. 

The newly developed General Amber Force Field (GAFF) was used to approximate the bonded interactions 
of all the simulation molecules, while the simulation package Spartan (Wave-function, Inc., Irvine, CA) was used at 
the Hartree-Fock/6-31g* level to obtain the necessary point charges for each of the atoms. To simulate a charging 
field, the charge on the carbon atoms of the graphite sheets was set to —0.0125 e per atom. The bonded (first three 
terms of Eq. ^) and non-bonded (last term) interactions in the AMBER Force Field are represented by the following 
potential-energy function: 
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where K^, Kg and Vn are the bond stretching, bending and torsional constants respectively, the constants A and B 
define van der Waals' interactions between unbonded atoms, and e is the electrostatic permittivity. The simulation 
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FIG. 2: Root-mean-square displacement of Lithium ions as a function of time. Diffusion is much faster with the additional 
oscillating electric field (amplitude 5 kCal/mol, frequency 25 GHz). 

package NAMD was used for the MD simulations, while the graphics package VMD [ll| was used for visualization 
and analysis of the simulation results. 

The model system representing the anode half-cell is composed of four graphite sheets (anode) containing 160 
carbon atoms each, two PF^ ions, and ten Li+ ions, solvated in an electrolyte made of 69 propylene carbonate and 87 
ethylene carbonate molecules (see Fig. [ija)). The graphite sheets were fixed from one side by keeping the positions 
of the edge carbon atoms fixed. 

B. Simulations and Results 

After energy minimization, the simulations were run at constant pressure using a Langevin piston Nose-Hoover 
method [13, Il3l | as implemented in the NAMD software package until the system has reached its equilibrium volume 
at a pressure of 1 atm and 300 K in the NPT (constant particle number, pressure, and temperature) ensemble. The 
system's behavior was then simulated for 200 ns (100 million steps) in the NVT (constant particle number, volume, 
and temperature) ensemble. Two observations were made: first, the Li"'' ions stayed randomly distributed within the 
electrolyte, and second, none of the Li+ ions had intercalated between the graphite sheets after 200 ns (see Fig.fTfb)). 

While the Lithium ions do not intercalate within the simulation time given above, it is expected that given enough 
time they will move towards the graphite sheets and get intercalated. To test whether intercalation is possible in such 
a model system, one of the Lithium ions was positioned between the graphite sheets at the beginning of a simulation, 
and we observed whether it diffused out from between the sheets. The Lithium ion stayed intercalated, even after 
400 ns. 

In order for intercalation to occur, the Lithium ion has first to diffuse within the electrolyte until it reaches the 
graphite electrode. Consequently, faster diffusion would result in faster intercalation and shorter charging time. In 
order to increase the diffusion of Lithium ions in the electrolyte, we explored a new charging method. In addition to 
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the charging field due to the fixed charge on the graphite carbons, an external oscillating square- wave field (amplitude 
5 kCal/mol, frequency 25 GHz) was applied in the direction perpendicular to the plane of the graphite sheets. Not 
only does this additional field increase diffusion, but also some of the Lithium ions intercalate into the graphite sheets 
within an average time of about 50 ns. Figure [2] shows a plot of the root-mean-square displacement of Lithium ions as 
a function of time for a system with and without an applied external field. The increased diffusion and intercalation 
indicate that a charging protocol involving an oscillating field may decrease the charging time and possibly increase 
the battery's power density. 



III. LATTICE-GAS MODELS OF CHEMISORBED SYSTEMS 



As mentioned in the Introduction, even the simplest electrosorption systems are extremely complicated. This 
complexity means that a comprehensive theoretical description that enables predictions for phenomena on macroscopic 
scales of time and space is still generally impossible with present-day methods and technology. (Note that MD 
simulations, such as those presented in SecHIl are only possible up to times of a few hundred nanoseconds.) Therefore, 
it is necessary to use a variety of analytical and computational methods and to study various simplified models of 
the solid-liquid interface. One such class of simplified models are Lattice-gas (LG) models, in which chemisorbed 
particles (solutes or solvents) can only be located at specific adsorption sites, commensurate with the substrate's 
crystal structure. This can often be a very good approximation, as for instance for halides on the (100) surface of Ag, 
for which it can be shown that the adsorbates spend the vast majority of their time near the four-fold hollow surface 
sites 114. A lattice-gas approximation to such a continuum model, appropriate for chemisorption of small molecules 
or ions jlSl. [l6l. il7.. il8>. .IS. i2Qi]. is defined by the discrete, effective grand-canonical Hamiltonian, 

(") 

HLG = ^h<i>^"^5]QC,]+H3-/i^C, . (3) 
n {ij) i 

Here, the lattice sites i are the preferred adsorption sites (the minima of the continuous corrugation potential), and 
Ci is a local occupation variable, with 1 corresponding to an adsorbed particle and to a solvated site. The sums 
Y^^/^j'^ and J2i '^un over all nth-neighbor pairs and over all adsorption sites, respectively, is the effective nth- 
neig hbor pair interaction, and runs over the interaction ranges. The term Ti.^ contains multi-particle interactions 
[2lL 122. l23l|. The sign convention is such that < implies repulsion, and jl > favors adsorption. Equation ((3|) is 
also easily generalized to multiple species (23 . [25l | . 

To connect the electrochemical potentials to the concentrations in bulk solution of species X, [X], and the electrode 
potential, E, one has (in the dilute-solution approximation) 

M,(r, [X],E) = + fcBTln([X]/[X]") - e r j^{E')dE' , (4) 

Jeo 

where is Boltzmann's constant, T the temperature, e the elementary charge, and jx{E) the electrosorption valency 
[H, [13, m, of X. The importance of the integral over the potential-dependent electrosorption valency (rather 
than just the product ej^(E)E) analogous to the case of potential-independent 7x) was pointed out in Ref. [30|. The 
quantities superscripted "0" are reference values that include local binding energies. The interaction constants and 
electrosorption valencies are effective parameters influenced by several physical effects, including electronic structure 



|2lL |22|. |23||. surface deformation, (screened) electrostatic interactions [31|, [32, |33| , and the fluid electrolyte |34l. |3 



The density conjugate to is the coverage relative to the number TV of adsorption sites, 

ex = A^"'^c, . (5) 



IV. CALCULATION OF LATTICE-GAS PARAMETERS BY DENSITY FUNCTIONAL THEORY 

There are many methods to estimate lattice-gas parameters. One of these is comparison of MC simulations (see 
Sec. IV]) of a LG model with experimental adsorption isotherms. For detailed descriptions of this method we refer 
to Refs. jlO, HH m, [13, [H, [39|. Here we instead concentrate on the purely theoretical method based on quantum- 
mechanical DFT calculations [23| . 

DFT is the most widely used method to calculate ground-state properties of many-electron systems. It is based 
on the Hohenberg-Kohn theorem, which states that all properties of the many-particle ground state can be expressed 
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(A) 



FIG. 3: (A) Cross section of a 3 x 3 supercell with O = 1/9. (B) Three-dimensional representation of the same cell and 
coverage. (C) Top view of a 3 x 3 surface and a 2 x 2 surface with various coverages. 



in terms of the ground-state electron charge-density distribution 40] and leads to the Kohn-Sham equations for 
single-particle wave functions [4T|. These are second-order differential equations, which include potential terms due 
to the ions and the classical Coulomb repulsive energy between the electrons, as well as the electronic exchange- 
correlation energy, and they are solved self-consistently. For surface structural studies, DFT is usually performed 
using pseudopotentials with slab models and plane-wave basis sets. The slab consists of a finite number of atomic 
layers, periodic in the direction parallel to the surface, which can either be repeated periodically in the third direction 
(separated by a vacuum interval), or not. The fluid solvent can be considered either as an effective continuum, or by 
molecular models. 

Here we present preliminary results on a DFT calculation of lateral interaction constants pertaining to a lattice- 
gas model for the adsorption of Br on single-crystal Ag(lOO) surfaces [1^, [13, [H, [s^ [13] ■ The lattice-gas model is 
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represented by Eq. ([3]) on a square lattice with lattice constant a = 2.95 A, Tis = 0, infinitely repulsive interactions 
for adparticles at nearest-neighbor sites, and the long-range repulsion 

(kij = 3 (t>nnn for nj > , (6) 

which is compatible with dipole-dipole interactions or elastically mediated interactions. (Here, rij is given in units of 
a.) Since the DFT calculations are performed in the canonical ensemble (fixed adsorbate coverage), p, in Eq. ([3]) is 
replaced by the binding energy of a single adparticle, i?b- 

We prepared slabs with seven metal layers, which were placed inside a supercell with periodic boundary conditions. 
Two different sizes of supercells were used: a 2 x 2 supercell with the size of 2a X 2a X 36.95 A, and a 3 X 3 supercell 
with the size of 3a X 3a X 36.95 A. The vacuum region above the surface was twice the thickness of the slab, and 
the orientation of the surface normal was in the z direction. One, two, and three Br atoms were placed on the 3x3 
surface to represent coverages & = 1/9, 2/9, and 1/3. Two Br atoms were placed on the 2x2 surface to represent 
O = 1/2, and one to represent = 1/4. Supercells with different coverages of Br are shown in Fig.lSl 

The DFT calculations were performed using the Vienna Ab Initio Simulation Package fVAS P') Esl . 0, \^ . The 
basis set was plane- wave, with the generalized gradient-corrected exchange-correlation function [i^. |47|| . and Vanderbilt 
pseudopotentials [i^ . The fc-point mesh was generated using the Monkhorst method [i^ with a 5 x 5 x 1 grid for the 
3x3 cells and a 7 x 7 x 1 grid for the 2x2 cells. All calculations were done on a 54 x 54 x 192 real-space grid. 

Individual DFT calculations provide total energies, E, and charge densities, p{x). The adsorption energy £^ads for 
a single adatom and the corresponding charge-transfer function Ap(x) are obtained from calculations of the adsorbed 
system and isolated slab and atoms as follows: 

Eads — [-Esyst — -Esiab] / A^'ads — Ebt (7) 

and [13] 

Ap{x) = [p(f)syst - P(^)slab]/Afads " p(^)Br , (8) 

where A^ads = A^0 is the number of adsorbed Br atoms in the cell, and the quantities subscripted Br refer to a single, 
isolated Br atom. 

Since the system is electrically neutral, the integral over space of Ap(a;) vanishes. The surface dipole moment is 
defined as 



p = y zAp{z)dz . (9) 



Kohn and Lau [51] have shown that the non-oscillatory part of the dipole-dipole interaction energy between adsorbates 
separated by a distance R behaves as 

2papt 

for large R (in our case larger than the nearest-neighbor distance). This result is twice what one might naively expect. 
Thus, the next-nearest-neighbor interaction constant from Eq. ^ would be 

2p^ 

(/-dip-dip nnn - 4^,^^3^^^ i^) 

with p obtained from the DFT by Eq. This estimate, which depends on Q, is included in Fig.|4]as solid circles. 

Alternatively, the interaction constant (/)nnn in the LG Hamiltonian, Eq. can be estimated by performing a 
nonlinear least-squares fit of the 0-dependent DFT adsorption energy i^ads hi Eq. ([7]) to 

Eads = — (/"nnn^e — -Eb© (12) 

with (^nnn — A{1 + BG)^, using the three fitting parameters A, B, and E\^. This is consistent with the theoretical 
prediction of Eq. pip with a dipole moment that depends linearly on O. The quantity 

Kj y 
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FIG. 4: Three different estimates of the LG interaction constant (pnnn- Circles: Based on Eq. (|10p with the dipole moment p 
directly obtained from the DFT calculation. Squares: Based on a three-parameter fit to the DFT adsorption energy i?ads as 
described in Eq. (|12|l . Diamonds: Based on minimizing mean-square deviations (MSD) from the estimate based on the DFT 
dipole moment p, constrained to retain a low value of the from the fit to -Bads. See discussion in the text. 



can be calculated numerically to any given accuracy for a particular coverage and adsorbate configuration. This 
estimate for 0nnn is included in Fig. 2] as solid squares. It does not agree particularly closely with the result obtained 
from the dipole moments. However, we found that the of the fit, considered as a function of the fitting parameters, 
was characterized by an extremely wide and shallow basin surrounding its minimum. We therefore further minimized 
the mean-square deviation (MSD) between the values of ^nnn obtained from this fitting procedure and those obtained 
directly from Eq. (|lip with the DFT values for p within the three-dimensional parameter region for which the original 
was close to its minimum. This procedure gave significantly improved consistency between the two estimates for 
0nnn, without a significant increase in x^. The final result is shown as solid diamonds in Fig. 01 and the corresponding 
parameters are listed in Table HI 

The average value of (/'nnn obtained by this method is consistent with that found by fitting equilibrium MC simu- 
lations (see Sec. IV A)) to experimental adsorption isotherms in aqueous solution (approximately —21 meV). However, 
no significant coverage dependence was found in the analysis of the experimental data [sol . [33 | . It is not surprising 
that results from in situ experiments and in vacuo DFT calculations should show some differences, and we find it 
encouraging that the average results are consistent. Application of the method described here to Cl/Ag(100) gave less 
consistent results than for Br, possibly indicating that the effective interactions for CI are not purely dipole-dipole in 
nature 15211. 
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TABLE I: Results for the fits of the B-dependent lattice-gas interaction constant ^nnn according to the two methods described 
in the text. Here, v is the number of degrees of freedom (number of data points minus number of parameters, here equal to 2) 
for the initial nonlinear least-squares fit of 0nnn to the DFT adsorption energy -Eads, while MSD is the mean-square deviation 
between this estimate and the estimate obtained directly from the DFT dipole moment. Minimizing MSD within the basin of 
low x'^ significantly reduces MSD (see the greatly improved agreement in Fig. |4ll without significantly increasing x^- 



Method 


A 


B 






MSD/u 


Min. 
Min. MSD 


-6.017 X 10*2 
-4.085 X IQ-^ 


-0.8632 
-0.7595 


3.102 
3.070 


2.362 X 10"^ 
2.675 X 10"^ 


1.803 X 10^* 
7.692 X 10"'^ 



V. MONTE CARLO SIMULATIONS 
A. Equilibrium Monte Carlo 

As a method to obtain equilibrium properties of a system described by a particular Hamiltonian, MC is more 
accurate than mean- field approximations, especially for low-dimensional systems near phase transitions (36l[53|. This 
is an effect of fluctuations which, while ignored or underestimated by mean-field methods, are very important in 
two-dimensional systems. Given the rapid evolution of computers and the relative ease of programming of MC codes, 
this is our method of choice for equilibrium and dynamic studies of both lattice-gas and continuum models. 

The goal of an equilibrium MC code is to bring the system to equilibrium as rapidly as possible, and then sample 
the equilibrium distribution as efficiently as possible. The only requirement is that the transition rates between two 
configurations c and c' satisfy detailed balance, 

n{c' c)/n{c -> c') = exp [- {H{c) - H{c')) /k^T] . (14) 

This result applies to both continuum and discrete systems, and Ti may be a classical potential of predetermined 
form, or the interaction energies can be calculated "on the fly" by DF T )54ll. The sampling can be accomplished 
with a number of different choices of the transition rates TZ{c' — > c) [H, l53l. l55l. [56l . ISTI . ISSI . ISQI . i&ol . [U, [63|, including 
Metropolis, Glauber, and heat-bath algorithms. It is important to note that the stochastic sequence of conflgurations 
generated by an equilibrium MC algorithm does not generally correspond to the actual dynamics of the system. 

B. Kinetic Monte Carlo 

To construct a MC algorithm producing a stochastic path through configuration space that is a good approximation 
to the actual time evolution of the system (in a coarse-grained sense), one can introduce transition states between 
the lattice-gas states. Only then can "Monte Carlo time," measured in MC steps per site (MCSS) in a lattice-gas 
simulation, be considered proportional to "physical time," measured in seconds (42| . In a Butler- Volmer approximation 
[H, HEl J the free energy of the transition state between lattice-gas configurations c and c' is given by 

n* (c, c') = A + (1 - a) Wlg (c) + an^G (c') , (15) 

where the symmetry constant a — 1/2 for diffusion but may be different for adsorption/de sorp tion 13611 . The "bare" 
barrier A must be determined by other methods. These may be ab initio calculations [H^ |64 IgsL [66| . MD 
simulations of the diffusion process on a short time scale as in Sec. in 3, i, i 0, i, or comparison of dynamic 
simulations with experiments [4^ . The most common choice of transition rate for KMC in chemical applications is 
the one-step algorithm [63, , 

7^ (c ^ c') = lyo exp [- {n*{c, c') - 7^lg(c)) /fceT] , (16) 

where is an attempt frequency (often of the order of a phonon frequency (10^ - 10^'^ Hz) , but see Ref. [43| for 
exceptions) that must be determined by other means. As we have shown previously [ssl. Is^ IstI. IssI. Is^. [60|. [61|. |62|| . in 
order to obtain reliable structural information from a KMC simulation, the transition rates must approximate the real 
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physical dynamics, which includes using transition states with proper energies. While the need for correct transition 
rates may seem obvious, it is regrettably often ignored in the literature. The most difficult barrier to estimate is that 
for adsorption/desorption, which requires reorganization of the adparticle's hydration shell. 

Since the transition rates used in KMC of activated processes are typically small, simulations that extend to 
macroscopic times must use a rejection-free algorithm, such as the n-fold way [69l F70| or one of its generalizations 
[63, Inl, [tS, [li m, [li, [li. Ill- These algorithms simulate the same Markov process as the "naive" MC approach 
of proposing and then accepting or rejecting individual moves. Although they require more bookkeeping (see the 
Appendix of Ref. [71} for an example), they avoid the large waste of computer time resulting from rejected moves. 



VI. ELECTROCHEMICAL FIRST-ORDER REVERSAL CURVE SIMULATIONS 



The First-order Reversal Curve (FORC) method was originally dev elop ed to enhance the amount of dynamic 
information extracted from magnetic hysteresis experiments [t^, [tI, [s^, l8l| . We recently proposed that the method 
can be further developed as an extension of traditional CV to study the dynamics of phase transitions in electrochemical 
adsorption [s^. [83|. 

This electrochemical FORC (EC-FORC) method consists of saturating the adsorbate coverage 8 in a strong positive 
electrochemical potential p, and, in each case starting from saturation, decreasing /I at a constant rate to a series of 
progressively more negative "reversal potentials" fir (see Fig.[5ja)). Subsequently, fl is increased back to the saturating 
p, at the same rate. (Saturation at negative potentials with reversal potentials in the positive range is also possible.) 
The method is thus a simple generalization of the standard CV method, in which the negative return potential 
is decreased for each cycle. This produces a family of FORCs, 0(/Ir,/Ji), where p,i is the instantaneous potential 
during the increase back toward saturation. In CV experiments, one actually records the corresponding family of 
voltammetric currents, 

■f- - dpi de {fir, Pi) ..„^ 

z(M.,M.) = -7e^^^, (17) 

where 7 is the electrosorption valency and e is the elementary charge (see Fig. [SJb)). 

The next step in extracting dynamical information from the FORCs or the corresponding currents is to calculate 
the FORC distribution, 

^_ 1 ^'^Q ^ 1 di{jir,Pi) 

2dprdpi 2je{dpi/dt) dpr 

This is shown in Fig. [SKc) as a contour plot commonly known as a FORC diagram in terms of the more convenient 
variables pb — {pr + Pi)/ 2 and pc — {pr ~ Mi)/2 [Z^, Geometrically, p is proportional to the vertical distance 
between adjacent current traces. 

To our knowledge, the data for our model of Br/Ag(100) [1^, [13, [H, js^ which are shown in Fig. O are 
the first FORC predictions for a continuous phase transition. All three panels are significantly different from the 
corresponding data for a discontinuous transition, such as seen in underpotential deposition (UPD). In particular, the 
FORC distribution for a discontinuous transition contains a negative region, while this does not appear for continuous 
transitions. (See details in [s^ . Issj.) Closely related to this negative region is an extremum of the current density 
during the return scan [sj]. EC-FORC analysis should be a useful and valuable method to distinguish between 
continuous and discontinuous phase transitions in experiments. 



VII. CONCLUSION 



In this chapter we have presented some applications of the statistical-mechanics based computer-simulation methods 
of Molecular Dynamics and equilibrium and kinetic Monte Carlo simulations complemented by quantum-mechanical 
density functional theory calculations of interaction energies. These include both highly technologically-oriented ap- 
plications to Lithium-battery technology, and basic-science investigations into adsorption on single-crystal electrodes. 
Our hope is that these examples and the list of references will encourage other workers in surface electrochemistry 
to take advantage of the recent spectacular advances in computational power and algorithmic sophistication to study 
ever-more detailed and accurate models of processes at solid-liquid interfaces. 
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FIG. 5: (a) Family of FORCs for our model of Br/Ag(100), corresponding to potential sweeps back and forth across the 
continuous phase transition between the disordered and c(2 x 2) phases. The bold arrows show the directions of the potential 
sweeps, and the vertical arrow indicates p,r for one of the FORCs. The bold curve is the FORC whose minimum lies closest to the 
critical coverage (shown in more detail in the inset). The thin curve in the middle is the equilibrium isotherm, (b) Voltammetric 
currents corresponding to the FORCs in (a), (c) Contour plot of the FORC distribution p, corresponding to the FORCs in (a). 
The jagged curve of dots in the upper part of the diagram corresponds to the minima of the positive-going curves in (a). The 
area above the curve corresponds to desorption, and the area below it to adsorption. The slanted, straight line corresponds to 
the bold curve in (a). After Ref. 
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